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ANALYSIS OF TURBULENT LIQUID-METAL HEAT TRANSFER IN CHANNELS 
WITH HEAT SOURCES IN THE FLUID - POWER-LAW VELOCITY PROFILE 

by Robert M. Inman 
Lewis Research Center 

SUMMARY 

An analysis is made to determine the heat-transfer characteristics for turbulent flow 
of a heat- generating liquid metal between parallel plates with wall heat transfer. The in- 
ternal heat generation is uniform over the channel cross section and along its length. 

The wall heat transfer is also uniform along the channel length. The analysis applies in 
the thermal entrance region of the channel as well as far downstream. The fluid is as- 
sumed to have a fully developed, turbulent power-law velocity profile which is unchanging 
throughout the length of the channel. The idealized eddy diffusivity profile proposed by 
Poppendiek is used. 

The solutions depend on the power -velocity exponent m and diffusivity parameter 

k = V 1 + 0. 01 i^PrRe^' where Pr is Prandtl number, Re is Reynolds number, and 
i ]s is defined as the average effective value of the ratio of the eddy diffusivity of heat 
transfer to that for momentum transfer. Numerical results for the wall temperature dis- 
tribution and Nusselt number variation are presented in graphical form for a 1/7-power- 
law velocity profile and for values of k ranging from 1 to 5. 

Results for the fully developed Nusselt numbers for liquid metal flow without internal 
heat generation are compared with existing calculations (based on the assumption that the 
eddy diffusivities of heat and momentum are equal) and exhibit good agreement. 


INTRODUCTION 

The study of flowing liquid metals with volumetric internal heat generation is of cur- 
rent interest in several sectors of modern technology. Such flows may occur, for ex- 
ample, in liquid metal magnetohydrodynamic generators for the generation of electric 
power in space and in liquid-metal-fueled nuclear reactors. In such devices, the liquid 
metal will be internally heated by the flowing electric current, by radioactive fission 


products, and, perhaps, by viscous dissipation. The fluid flow in these devices may, in 
addition, be accompanied by wall heat transfer. 

This investigation is concerned with fully developed turbulent channel flows in which 
the internal heat generation per unit volume of liquid is uniform. A parallel- plate chan- 
nel, which approximates a large-aspect-ratio channel, is analyzed in this work as it has 
applications in the aforementioned devices. 

Attention is focused here on the case where the wall heat flux is uniform along the 
length of the channel. A factor of importance for the proper operation of these devices is 
that of maintaining a satisfactory temperature distribution along the channel walls. The 
designer, therefore, must be able to compute the wall temperature under the conditions 
of internal heat generation and wall heat transfer. The problem involves studying the ef- 
fect of internal heat generation and wall heat transfer on the wall temperature distribu- 
tion. 

The heat-transfer behavior for this type of situation has been considered in two pre- 
vious papers. In reference 1, the effect of an internal heat source on heat transfer in 
round tubes and flat ducts is considered under the conditions of a uniform velocity profile 
(or slug flow) and turbulent heat transfer occurring solely by molecular conduction. Re- 
sults were obtained in both the thermal entrance and fully developed regions. 

The solution in reference 2 is for liquid-metal flow in a parallel- plate channel. The 
limitation to heat transfer by molecular conduction only is eliminated, and the effect of 
transverse thermal eddy diffusivity variations is included through consideration of a sim- 
plified eddy diffusivity profile. Results were obtained in both the thermal entrance and 
fully developed regions. The established turbulent velocity profile was again approxi- 
mated by a slug-flow profile. 

In the present analysis, this last limitation of uniform velocity profile, or slug flow, 
is also eliminated. The established turbulent velocity profile is represented by a power- 
law expression. Together with a parameter (to be presented later) that characterizes the 
effect of eddy transfer at moderate Prandtl or high Reynolds numbers, the wall tempera- 
ture distribution and heat-transfer characteristics can be obtained for any value of the 
exponent in the power-law velocity expression. A value of the exponent of 1/7 is used in 
the evaluation of the solutions. 


ANALYSIS 

A necessary prerequisite for the solution of the turbulent convective heat-transfer 
problem is a specification of the velocity profiles and the heat transfer eddy diffusivity 
profiles. Attention, therefore, will first be directed to the velocity and eddy diffusivity 
problems. 


2 


Velocity Distribution 


The simplification provided by the slug-flow assumption in references 1 and 2 made 
it possible to obtain exact mathematical solutions to the governing energy equation. The 
actual turbulent velocity profile in channel flow is not uniform but, instead, falls off to- 
ward zero near the walls. In an attempt to simplify analyses and, yet, deal realistically 
with turbulent flows confined in pipes and channels, much attention has been given to 
power laws (ref. 3). In order to approximate in a simple manner the velocity distribution 
for turbulent viscous flow in a pipe, a power-law velocity distribution of the type 



is frequently considered. The exponent m in this power-law velocity expression has 
been experimentally determined to be 1/7 to a Reynolds number of approximately 10^. 

For Reynolds numbers larger than 10^, a better approximation is obtained by the power 
one- eighth, one-ninth, or even one-tenth. Therefore, the slug-flow velocity profile is 
expected to be a good approximation for the actual turbulent flow at very high Reynolds 
numbers for which the velocity distribution is extremely flat. 

Although the power-law velocity expression was developed for flow in circular pipes, 
subsequent study revealed that it represents satisfactorily the established turbulent veloc- 
ity distribution for flow between parallel plates (ref. 3) when written in the form (fig. 1) 

u = A(l-X) m = A(l-y.) m (2) 

The mean fluid velocity u of the parallel-plates system can be expressed in terms of the 



Figure 1. - Geometry and coordinate system. Turbulent 
velocity profile represented by power-law distribution. 
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hydrodynamic parameter A as follows: The mean velocity is given by the definition 



Substituting equation (2) into equation (3) and carrying out the integration yield 

A = u(m + 1) 

Therefore, equation (2) can be rewritten as 

U(y’) = = (m + 1)(1 - y') m 

u 

With an exponent m of 1/7, equation (4) yields 

U(y') =- (1 - y’) 1/? 

7 


while an exponent of zero in equation (4) yields the uniform velocity profile 


( 3 ) 


(4) 


(5) 


U(y’) = 1 (6) 

Both the 1/7-power-law velocity and the uniform-velocity profiles have been used to gain 
an understanding of turbulent heat-transfer characteristics for liquid-metal flow in tubes 
and channels in the absence of internal heat generation (refs. 4 to 11). Both profiles have, 



Dimensionless radial 
distance, y/a 

Figure 2. - Velocity profile 
for fully developed turbu- 
lent flow between parallel 
plates. 
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in addition, been used in the analyses of hydromagnetic effects in turbulent channel flows 
(refs. 12 to 15). The assumption of a uniform velocity distribution, for example, is be- 
lieved reasonably satisfied for a large magnetic-field strength; in practical design calcu- 
lations, the nonuniformity of the velocity profile would need to be considered when esti- 
mating the heating of the fluid caused by ohmic and viscous losses. The velocity profiles 
represented by equations (5) and (6) are shown in figure 2. 


Eddy Diffusivity Distribution 

In most of the previously mentioned analytical studies of turbulent liquid-metal heat 
transfer, it is postulated that the thermal eddy diffusivity is small compared with the 
thermal molecular diffusivity and may be neglected; therefore, radial or transverse heat 
transfer is by molecular conduction only. This is believed a reasonable postulate when 
the product of Reynolds and Prandtl numbers is below a specified limit (about 100). The 
practical design of moderate Prandtl or high Reynolds number liquid-metal systems, 
however, where the effect of radial or transverse eddy diffusivity variations must be in- 
cluded, has prompted the development of idealized eddy diffusivity functions. Poppendiek 
(refs. 9 and 10) has analyzed turbulent heat transfer for liquid-metal flow between paral- 
lel plates for the uniform wall-temperature case and with no internal heat generation by 
considering a simplified eddy diffusivity function. This function was also used in refer- 
ence 2. The solution is shown in references 2, 9, and 10 to reduce correctly to known 
specific solutions of the general case. 

The idealized eddy diffusivity profile proposed by Poppendiek (refs. 9 and 10) is used 
in the present analysis to account for transverse heat flow consisting of eddy transfer as 
well as molecular conduction. This simplified profile approximates the actual one in the 

O e 

regions nearest the wails over a Reynolds number range of approximately 5x10 to 1x10 . 
For symmetrical wall heating conditions considered herein, the use of this idealization 
should introduce little error. 

The idealized momentum eddy diffusivity distribution varies linearly with distance 
from the channel centerline and as the nine-tenths power of the Reynolds number: 

— = 0. 01 Re 0 - 9 (1 - y') (7) 

v 

In anticipation of a later need, it will be convenient to introduce the thermal eddy diffu- 
sivity parameter Xrp defined by 
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(8a) 


\ T = Pr — = 0. 01 i^PrRe 9 ' 9 (1 - y') 

v 

where is the ratio of the eddy diffusivity of heat to that of momentum. No generally 
accepted relation between the two quantities has yet been established. It has been found 
convenient, however, to consider the parameter , which is the effective average value 
of across the channel (ref. 16). The parameter X T is then given by 

X T = 0. 01 ^PrRe°* 9 (1 - y’) (8b) 

It is convenient to write the parameter X,p in terms of a diffusivity parameter k defined 
by 


k 2 h 1 + 0. 01 ^PrRe°‘ 9 (9) 

The parameter X,p then takes the form 

X T = (k 2 - 1)(1 - r) (10) 

It is noted that heat transport in the fluid occurs solely by molecular conduction (i. e. , 

Xrp = 0) for y* = 1 (i. e. , at the channel walls) or for k = 1. It is clear from equation (9) 
that k can never physically assume a value less than one. 

Now that the velocity and eddy diffusivity distributions have been specified, the solu- 
tion of the heat-transfer problem is undertaken. 


Heat -Transfer Problem 

The geometry of the parallel -plate channel and coordinates are shown in figure 1. 
The fluid possesses a fully developed velocity profile which is unchanging with length. 
The fluid enters the channel at x = 0 with a uniform temperature at the value t^ and is 
heated both by the internal heat generation in the fluid and by the heat flux at the channel 
walls. An analysis is performed to determine the resulting wall temperature and heat- 
transfer characteristics as functions of the axial position. The liquid metal is assumed 
to have constant physical properties, and only steady- state flow and heat transfer are 
considered. 

The energy equation for the fluid temperature is 
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Q 


( 11 ) 


u 


9t 

3x 


d_ 

9y 


(« + e H ) — 
9y 


Pc T 


where a and e„ represent, respectively, the molecular and eddy diffusivities for heat. 
To obtain the energy equation in this form, it has been assumed that viscous dissipation 
and axial heat conduction are both negligible compared with heat conduction in the trans- 
verse direction. This assumption has been shown (ref. 17) to introduce a negligible error 
for Pe > 100. 

It is convenient to write the temperature t as the sum of the temperature for the 
problem where there is a uniform volumetric internal heat generation Q in an insulated 
channel tQ, and the temperature for the problem where there is a wall heat transfer q w 
without internal heat generation t^. Then the temperature in the combined problem is 
simply given by 


‘ - *Q + ‘q 


( 12 ) 


The governing equations and boundary conditions for tQ and t may be written as 


and 


at. 


Q _ 9 


9x 9y 


(a + e H ) 


at. 


Q 


9y 


_Q_ 

Pc„ 


(13a) 


9t, 




9y 


^ = 0 at y = a (insulated wall) 


9t, 


9y 


^ = 0 at y = 0 (symmetry) 


t„ = 0 at x = 0 (entrance condition) 

^ j 


(13b) 


u 


dx 0y 


at a 

(a + e H ) — S 


9y 


(14a) 
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3t q 

— - = — at y = a (specified heat flux) 

ay k 


at 


ay 


2=0 at y = 0 (symmetry) 


tq = tj at x = 0 (entrance condition) 


(14b) 


Before attempting to solve for the temperatures tQ and t , it is convenient to write the 
preceding equations as follows: 


and 


at, 


u 


Q _ a 


ax’ ay' 


at, 


Q 


(i + x T ) — 
ay’ 


Qa* 


(15a) 


at, 




ay' 


Q_ 


= 0 at y' = 0 and y' = 1 


tQ = 0 at x' = 0 


(15b) 


at 


u 


q _ 3 


ax* ay' 


at' 

(i + x ) -2 
1 ay 


(16a) 


at 


q _ ^w £ 




ay' k 


at 


at y' = 1 


— 2 = o at y’ = 0 ^ 

ay 


t = t, at x’ = 0 

4 •*- 


j 


(16b) 


Equation (10), for X T , may be inserted into equations (15a) and (16a). It is also conven- 
ient to introduce the change of variable 


which then allows equations (15) and (16) to be written more compactly as 


and 


v _ a o L (r^Q) , Qa2 

9x’ r 9r y 9r J k 

dt O 

— — =0 at r = 1 and r = k 

9r 

tq = 0 at x’ = 0 


(18a) 


(18b) 



_9_ 

9r 



(19a) 



9r 


1 V 



at r = ll 


—2=0 at r = k 
9r 


t 


q 



at x' = 0 J 


(19b) 


2 * 

where a Q = (k - 1) /4. The power-law velocity expression, in terms of the new variable 
r, is given by 


U(r) 


m + 1 ( r 2 

o m 

(k 2 - 1) 


m 

1) 


( 20 ) 


Equations (18) and (19) complete the formulation of the boundary- value problems. 

The problem of solving for t^ and t q is most conveniently attacked by separation of en- 
trance region and fully developed solutions: 



I 


t - t , + t 
q q, d + q, e 


(21b) 


The two problems as defined by equations (21a) and (21b) are analyzed separately. The 
solutions are then combined to yield the results for the general case of both internal heat 
generation and wall heat transfer. 

The solution for t~. which applies both in the entrance and fully developed regions is 



K 


- x r + F(r) + 



A n R n (r) e 



( 22 ) 


in which the radial function F(r) is the solution of the mathematical system 

[U(r) - 1] (23a) 

r dr \ dr/ a Q 


0 at r = 1 and r = k (23b) 

dr 



U(r)F(r)r dr = 0 


(23c) 


The condition at r = k is automatically satisfied when the condition at r = 1 is satisfied. 

The last condition is the result of the consideration of an overall energy balance on the 

o 

fluid in the fully developed region. The a* and R n are, respectively, the eigenvalues 
and eigenfunctions of the Sturm- Liouville mathematical system 


1 d_ 
r dr 



+ — U(r)A = o 
a o 

> 


dR n 


= 0 


dr 


at r = 1 and r = k J 


(24) 


The coefficients A are determined to satisfy the condition that t„ = 0 at the channel 
n W 

entrance: 
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( 25 ) 


A n 


/* k 

U(r)F(r)R n (r)r dr 
U(r)R^(r)r dr 


Equation (23a) for the radial function F(r) can be integrated directly. The resulting 
expression for F(r) is 


where 


F(r) = — 
a„ 


I (r) - < r : .i) + A In r 
1 4 2 


+ F(l) 


(26) 


F(l) = 


1 

~ I 2 (k) - I 3 (k) - | I 4 (k) 

(27a) 

O 

rt 

i 


/» r /•? 



i f U(o)o' do 


= ) 

r 1 d/3 (27b) 

P 




= j 

/* k (r 3 - r)U(r)dr 
1 

(27c) 

i 3 « = _ 

/* k 

J U^I^rJr dr 

(27d) 

i 4 « = 

/* k 

/ U(r)r In r dr 

(27e) 


and o' and /3 are dummy variables. 

Similarly, the solution for t is 

VI 
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( 28 ) 


t " 

-3 - = x' + G(r) + 

3w a 


OO 

I - 


~ a t x ' 

( r )e n 


n=l 


in which the radial function G(r) is the solution of the mathematical system 


1 A ( r ) 

r dr \ dr/ 

1 = — U(r) 
a o 

(29a) 

dG = _ 1 

at r = 1 

(29b) 

dr A 

dG=0 

at r = k 

(29c) 

dr 

/* k 

/ U(r)G(r)r dr = 0 

(29d) 


As in the foregoing, the a ^ and R n (r) are the eigenvalues and eigenfunctions of equa- 


tion (24), but now the coefficients B n are 


B n=“ 


/* k 

U(r)G(r)R n (r)r dr 
U(r)R^(r)r dr 


(30) 


in order to satisfy the condition t^ = tj at the channel entrance. 

Equation (29a) for the radial function G(r) can also be integrated directly to yield 


G(r)=i 


a o L 


Ij(r) - In r 


+ G(l) 


(31) 


where 
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■ in 


ii i 


i ii 


i i i 


it r 



G(l) = - 


1 


(32) 


! 3 ( k ) - I 4 (k)] 

To complete the solution, it is necessary to find the eigenvalues and eigenfunctions 
of equation (24) for specified values of the exponent m in the power- law velocity expres- 
sion, and to carry out the integrations for the A R and B R in equations (25) and (30). 

The determination of these quantities will be considered later. 

Combining the results of the previous paragraphs, it is possible to write the solution 
for the situation where internal heat generation and wall heat transfer occur simultane- 
ously. The temperature distribution for this combined problem, according to equa- 
tion (12), is found by adding the contributions due to each of the separate problems. The 
temperature difference t(x',r) - t^, which applies both in the entrance and fully developed 
regions, therefore, can be written as 



t-V 


. Qa 


Z -a x T q a \ ^ -crx f 

A n R n^ e n x’ +G(r)+ \ B n R n (r)e n 

n=l n=l 


(33a) 


It is convenient to rewrite equation (33a) in the nondimensional form 


t - t. 


- = (1 + R)x' + F(r) + 


Qa 


w o 

Z -ax 1 

A nV r > e n +R 


n=l 


G(r) + 


9 

B n R n (r)e n 


n=l 


(33b) 


where 


R = 


u w 

Qa 


(34) 


The heat-flux parameter R is essentially the ratio of the heat transferred at the channel 
walls to the heat generated internally. 

Of particular practical interest is the wall temperature variation corresponding to 
internal heat generation and wall heat transfer. This quantity can be found from equa- 
tion (33b) by evaluating the equation at r = 1: 
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(35) 


t - t. 

— = (1 + R)x T + F(l) + 

Qa 2 


OO o °Q 

Z — 

A n e n + R G(l) + \ 


n=l 


B n e 


-a 2 x'| 
n 


n=l 


Another result of practical engineering importance is the variation in the wall- to bulk- 
temperature difference. The bulk temperature t^ is given by 


t b (x') = tj + 


I Qa“ + V' 


(36a) 


or, alternately, 


t - t. 

— = (1 + R)x’ (36b) 

Qa 2 

K 


The bulk temperature rises in a linear fashion along the channel. The difference between 
the wall and bulk temperatures at all stations along the channel is given by 


- y . ' .- 1 ? = f(i) + 

Qa 2 



A n e 


a l x ' 


+ R 


G(l) + 



(37) 


A useful form of the wall temperature results, given by equation (37), is obtained by sep- 
arate consideration of the insulated or adiabatic wall (R = 0) and nonadiabatic wall (R * 0) 
cases. 

The wall- to bulk-temperature difference for the insulated wall case, denoted by 

(t - t b ) , is found by evaluating equation (37) with R set equal to zero: 
w D q 


(t W - 
Qa 2 


= F(l) + 



(38) 
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By combining equations (37) and (38), the following expression for the ratio (t - t^)/ 
(t w - t b ) is obtained: 



n=l 


The ratio given by equation (39) depends separately on the velocity distribution and the 
diffusivity parameter, as well as on x\ The departure of (t w - t^)/(t - t^) from unity 

is a measure of wall heat-transfer effects. 

The wall- to bulk-temperature difference for the fully developed situation (x' — °°) is 
found from equation (37) to be 

(t w " V . 

= F(l) + RG(1) (40) 

Qa 2 

K 

A convenient rephrasing of equation (37) may then be carried out by introducing the fully 
developed wall- to bulk- temperature difference, which yields 



The departure of (t - t.) / (t - t, ) from unity is a measure of the thermal entrance 

W D W D ^ 

effects. 

The heat-transfer characteristics of both laminar and turbulent flow are customarily 
represented in terms of a heat-transfer coefficient h s q^^/(t - t^) and a Nusselt num- 

ber Nu s hDjj//c, where Djj is the hydraulic diameter (Djj = 4a for the parallel-plate 
channel). With these definitions, it follows from equation (37) that 
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I 



Nu = 


4 


(42) 


where 


r 


F(l) + 



_ -Of 

V 


n 


t 


n=l 


+ G(l) + 



_ -Of 

B n e 



S 


Qa _ 1 

a R 
H w 


(43) 


The Nusselt number depends on the important parameters k and m, the longitudinal 
position x', and the inverse heat-flux ratio S. The Nusselt number in the absence of in- 
ternal heat generation, denoted by Nu Q , is obtained by evaluating equation (42) with 
S = 0. The result is given by 


Nu q = 


G(l) + 



(44) 


As a matter of general interest, the Nusselt number Nu Q in the fully developed region, 
denoted by Nu Q d , is considered. This is obtained from equation (44) by considering the 
limit x' — 


Nu 


0, d 


4 

G(l) 


(45) 


Use will be made of the fully developed Nusselt number in the section entitled RESULTS 
AND DISCUSSION. 


Limiting Case: Diffusivity Parameter k = 1 

The analysis of the preceding section led to the determination of the wall temperature 
distribution and heat-transfer characteristics with transverse heat diffusion occurring by 
turbulent eddying as well as by molecular conduction. It is of interest to examine the 
limiting case where the thermal eddy diffusivity is negligible compared with the thermal 
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molecular diffusivity. In this situation, transverse heat transfer is by molecular conduc- 
tion only; this corresponds to the limiting value for k of one. A number of equations 
presented in the preceding section assume indeterminate forms for the limiting case 
k = 1, and the evaluation of these results is a problem of considerable mathematical com- 
plexity. It is more convenient in this connection to reconsider equations (15) and (16) and 
to set Xrj, = 0 therein. The solution of the subsequent boundary value problems then 
yields the temperature distribution and heat-transfer characteristics in the absence of 
heat conduction by turbulent eddying. The methods of solution are similar to those out- 
lined in the earlier portion of this investigation; therefore the details of the derivation 
are omitted. The solutions for the wall- to bulk-temperature difference t - t^ and the 
bulk temperature change t^ - t. are 



Qa 2 


K 


F(l) + 



(46) 


t H - t, 

— = (1 + R)x' 

Qa 2 

K 


(47) 


in which 


F(l) = I 1 (l)-i[l-^(l)]-Ig(l) 
2 

(48) 

G(l) = I X (1) - Ig(l) 

(49) 

I x (l) = J* | ' ^ U(<*)da d| 

(50a) 

^(i) = jT u(m 2 d^ 

(50b) 

Ig(l)= / U^Ij^dT 

(50c) 
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and 


^(r) = 


■/[/ 


U^da? 




(50d) 


The variables a, %, and r are dummy variables. The coefficients C = C Y (1) and 
— „ „ ,,, n n n' ' 

D n = D n Y n (l) are given by 



(51) 



(52) 


o 

The /3 fi and Y n appearing in equations (46), (51), and (52) are, respectively, the eigen- 
values and eigenfunctions of the mathematical system 


d Y 

5 + U(y')A = 0 

9 n n 

dy’ 2 


^n 

dy’ 


= 0 at y' = 0 and y f = 1 

-J 


(53) 


The important wall temperature ratios (t w - t b )/(t w - t fe ) and (t w - t fe )/(t w - t b ) are 

0 u d 

readily obtained from equation (46) as 
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Ill 


G(l) + 


V^b 

■ V, 


oo o 


= 1 +R 


n=l 


QQ o 

Z 


(54) 


F(l) + 


n=l 


F(l) + 


V ~ 


<V-*b> 


-V=1 


QQ n 

Z v* 


n=l 


4 X ' 

n + R 

OO rt 

V s - -/3“x' 

G(l) + \ D n e n 


n=l 


(55) 


d F(l) + RG(1) 

The Nusselt number in the present situation is given by the following equation: 

4 


Nu = 


F(l) + 


n=l 


(56) 


oo 9 oo 9 

Z - "^n ~ - -0JV 

C n e n + G(l) + 2^ D n e n 


n=l 


Then the Nusselt number in the absence of internal heat generation is obtained from 

4 


Nu 0 = 


(57) 


G(l) + 


OO 9 

Z _ -rx' 


n=l 


from which it follows that 


Nu 


0, d 


4 

G(l) 


(58) 
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Transverse Distribution Functions R(r), Y(y') 


Attention is now directed to the Sturm- Liouville eigenvalue problems (eqs. (24) and 
(53)). Consideration will be given first to the determination of the function R(r) which is 
the solution to equation (24). To obtain a solution of equation (24), it is necessary that 
the variation of U with r be specified. When the velocity profile (eq. (20)) is intro- 
duced into equation (24), the governing equation and boundary conditions for R(r) are 


1 d_ 
r dr 



+ ^ r2 - D 


m 


R 


n 


0 


> 


dR 

— - = 0 at r = 1 and r = k 
dr 




(59) 


where 


x 2 _ 4 < m + 1 ) g n 

n 9 m+2 

Or - i) 


A different set of eigenvalues and eigenfunctions will be obtained for each independent 
value of k and m. For m = 0 there results a Bessel equation which, in general, has 
solutions in terms of ordinary Bessel functions of integral order. For this case, the 

n 

first five eigenvalues cT and coefficients B n have been given in reference 2. (For 

m = 0, the coefficients A n are equal to 0). For positive, fractional values of m, an 
analytical solution for equation (59) in terms of tabulated functions for arbitrary values of 
k does not appear possible. Therefore, it was decided to integrate equation (59) numer- 
ically for a value of m of 1/7 by the Runge-Kutta method and to determine the desired 
eigenvalues by trial and error. In addition, the coefficients A r and B n , given by 

1) 

R n (s)ds 

(60) 

9 ^9 /*(k “ 1) rn 2 

(k 2 - 1) A 2 f s m R 2 (s)ds 
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(61) 


B n = 


4r2( 0 ) 


<k‘ 


-<x 


(k 2 -l) 


m 


R 2 (s)tis 


2 

where the change of variable r - 1 = s has been introduced, are to be obtained by nu- 
merical integration. Solutions were carried out for values of the diffusivity parameter 
k of 1.25, 1.67, 2.50, and 5. 00. The first four eigenvalues ex 2 and coefficients A n 
and B n thus obtained are listed in tables I and n, respectively. 

TABLE I. - LISTING OF EIGENVALUES 
a 2 n FOR 1/7-POWER-LAW 
VELOCITY EXPRESSION 


Diffusivity 

Eigenvalue 

parameter, 

2 

2 

2 

2 

k 

a l 

a 2 

a 3 

« 4 

1.25 

13.06 

51.25 

114.2 

203. 1 

| 1.67 

18.89 

73.25 

163.0 

288.9 

2. 50 

34.95 

131. 5 

290.9 

514. 5 

5. 00 

115.8 

412.0 

894.9 

1566 


TABLE II. - LISTING OF COEFFICIENTS FOR 1/7-POWER- LAW VELOCITY EXPRESSION 


Diffusivity 




Coefficient 




parameter, 









k 

- A i 

“ A 2 

~ A 3 

- A 4 

- B i 

- B 2 

- B 3 

~ B 4 

1.25 

0.01379 

0. 001478 

0.0005431 

0. 0002160 

0. 18784 

0.05045 

0.02322 

0.01336 

1.67 

. 01052 

.001071 

. 0004058 

.0001571 

. 15035 

.04187 

.01945 

.01124 

2. 50 

.00659 

. 000646 

.0002561 

.0000956 

. 10072 

. 03053 

. 01427 

. 00848 

5.00 

. 00237 

. 000238 

.0001021 

.0000365 

. 04009 

.01484 

. 00764 

. 00463 


There remains to consider in this section the eigenvalues and eigenfunctions of equa- 
tion (53) and to carry out the evaluation for the coefficients C n and D n in equations (51) 
and (52). When the velocity profile (eq. (4)) is introduced into equation (53) and when the 
change in variable r] = 1 - y' is introduced, the governing equation and boundary condi- 
tions for Y (tj) are found to be 
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( 62 ) 


d 2 Y. 


d Irj 


n m, 2 V n 
— + 77 toY = 0 
i ' n n 






dY 

=0 at 7] = 0 and r] = 1 

dr? 


2 2 

where u> n = (m + 1)0^. A different set of eigenvalues and eigenfunctions will be obtained 
for each value of m. For the limiting case m = 0, a solution of equation (62) that satis- 
fies the boundary conditions is found readily to be 


Y n (r?) = E 1 cos [ n7r ( 1 “ V)] 


= Ej cos n7ry’ 


(63a) 


2 2 

where the eigenvalues w n = /3^ are given by 


= ( m ) 2 


(63b) 


The result is in agreement with that in reference 1. The constant Ej is an arbitrary 
constant. It is readily shown that, in addition, the coefficients D n are given as 


D n = - 


(63c) 


(™y 


The solution to equation (62) for positive fractional values of m remains to be deter- 
mined. It can be shown that the general solution to the differential equation (62) is 


Y(r?) = 7? 1 / 2 




E 0 J 


) + E 3 Y l/(m+2)(^7J 


(m+2)/2i 


(64) 


where J;iy( m+ 2) an d Y l/(m+2) are m + 2) order Bessel functions of the first and 

second kind, respectively, and Eg and Eg are arbitrary constants. The solution of 
equation (64) involves Bessel functions of fractional order. The eigenvalues a; 2 are ob- 
tained by requiring the general solution to satisfy the boundary conditions. If character- 
istic numbers and zeros were available for these fractional- order functions, then an an- 
alytical solution for the important eigenvalues and coefficients would be possible. In 
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particular, a value of the exponent m of 1/7 is to be used in the present evaluation of the 
solution. However, the author was unable to find any tabulation of Bessel functions (and 
derivatives) of the fractional order 7/15; an analytical solution for Yfa) was, therefore, 
not possible and a numerical technique (Runge-Kutta method) had to be employed. 

The first four eigenvalues /3^ = (7/8)u>^ and coefficients C ft and D n , alterna J ^ly 
given as 


C 


n 


D 


n 


V°> X' Y n <7,) *’ 


Y>) 






( 65 ) 


( 66 ) 


were evaluated by utilizing numerical integration. The eigenvalues and coefficients are 
listed in table HI. 


TABLE HI. - LISTING OF EIGENVALUES 
AND COEFFICIENTS FOR 1/7- POWER- 
LAW VELOCITY EXPRESSION 


Index, 

n 

Eigenvalue, 

/3 2 

'n 

Coefficient, 

-c n 

Coefficient, 

- 5 n 

1 

10.22 

0.0163500 

0.21510 

2 

40.31 

.0018452 

.05702 

3 

90.23 

. 0006639 

.02615 

4 

160.0 

.0002694 

.01503 


Constants F(l), G(l), F(l), and G(l) 

To complete the solution for the wall temperature distribution and heat-transfer 
characteristics it is necessary to carry out the integrations for F(l) and G(l) in equa- 
tions (27a) and (32) and for F(l) and G(l) in equations (48) and (49). Some of the integra- 
tion is elementary; the remainder requires changes of variable or repeated integration by 
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parts. The considerable amount of analytical detail is omitted, with the results pre- 
sented as follows: 


F(l) = 


m + 1 _ 4(m + 1) 


/ o\ai 2 o 2m+3 

(m + 2)(k - 1) ^2 _ 


/'V 


(T - 1) 

£ 


m+l 




r dr 


In k + 


(k 2 - 1) (k 2 - 1) 


-i- / 

m+3 / 

- 1) *1 


to - d* (67a) 

<P 


G(l) = 


2 Ink- 4 < m + 1 ) 

/i 2 1X o 2m+3 

(k - !) (k 2 - !) 


/H/ 

/ 

9 m+2 / 

nr - 11 •'l 


(r - 1) 


m+l 


d£ 


r dr 


(ff...-. 1 ) m+1 d<p (67b) 




F(l) = 


m + l m + l 


2(m + 3) (m + 2) (2m + 3) 


(67c) 


G(l) = 2 ( m + 1 ) 2 - (67d) 

(m + 2) (2m + 3) 

For the special case m = 0 (uniform velocity profile), it is easy to verify that these con- 
stants reduce to 


F(l) = F(l) = 0 


2k 2 In k - (k 2 - 1) - — (k 4 - 1) 


G(l) = 


(k 2 - 1)' 


(68a) 


(68b) 
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1 

G(l) = — 

3 


These results have been obtained earlier in reference 2. For positive fractional values 
of m, an analytical solution for F(l) and G(l) is not possible because the integrals con- 
tained therein are intractable, and numerical integration has to be employed. Solutions 
were carried out on an electronic computer for a value of m of 1/7 and for parametric 
values of k of 1.25, 1.67, 2.50, and 5. 00. Numerical values of F(l) and G(l), in ad- 

/n/ /v 

dition to F(l) and G(l) evaluated from equations (67c) and (67d) for m = 1/7, are listed 
in table IV. 


TABLE IV. - LISTING OF CONSTANTS 
FOR 1/7- POWER- LAW VELOCITY 
EXPRESSION 


Diffusivity parameter, 

Constant 

k 

F(l) 

8(1) 

1.00 

0.01950 

0. 3711 


F(l) 

G(l) 

1.25 

0.01632 

0. 3263 

1.67 

. 01236 

.2665 

2. 50 

. 00705 

. 1877 

5. 00 

. 00279 

.0867 


RESULTS AND DISCUSSION 

With the numerical information in tables I to IV, the longitudinal variations of the 
dimensionless adiabatic or insulated wall temperatures were evaluated from the analyti- 
cal solution (eq. (37)) and the results are plotted in figure 3. In interpreting this figure, 
it is important to note that the Reynolds and Prandtl numbers appear both in the param- 
eter k and in the abscissa. The information given in this plot permits evaluation of the 
insulated wall- to bulk-temperature difference at various stations along the channel as a 
function of Reynolds and Prandtl number. These results apply for a value of the power- 
law velocity exponent of 1/7. 

Another quantity which is of practical interest is the fluid bulk temperature variation 
along the length of the channel. The bulk temperature is given relative to the tempera- 
ture of the fluid at the entrance to the channel. The dimensionless bulk temperature 
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Dimensionless longitudinal distance, (x/a)/Pe 

Figure 3. - Wall temperature results for inter- 
nal heat generation in channel with insula- 
ted walls. Power-law velocity exponent, 1/7. 



Dimensionless longitudinal distance, (x/a)/Pe 

Figure 4. - Bulk temperature results for 
internal heat generation with wall heat 
transfer. 


given by equation (36b) has been plotted in figure 4 for parametric values of the heat-flux 
parameter R. 

Positive and negative values of the parameter R are considered in the figure. In 
the present analysis, Q is taken to be positive (a heat source). A positive value of R, 
therefore, implies that q w is positive, or that heat is being transferred from the walls 
to the fluid. A negative value of R, on the other hand, implies that q w is negative and, 
therefore, that heat is being transferred from the fluid to the walls. The case of R = 0 
corresponds to the thermally insulated wall. For the special value of R = -1. 0, the fluid 
flows with a bulk temperature which is unchanging with length. For any other arbitrary 
value of R, the bulk temperature rises (or falls) in a linear fashion along the channel. 

The ratio of nonadiabatic to adiabatic wall- to bulk-temperature differences (eq. (39)) 
is plotted in figure 5 as a function of the dimensionless axial distance along the channel 
for parametric values of R using the numerical data listed in tables I to IV. The infor- 
mation given on these plots, used in conjunction with figure 3, permits evaluation of the 
wall- to bulk-temperature difference at various stations along the heated channel walls. 

It is noted that, for some negative values of the heat flux ratio R, the wall- to bulk- 
temperature difference may be negative, which means that t fe is larger than t . This 
is understandable, however, if it is recalled that t is a local value along the wall, 
while t^ is an average value over the entire cross section. 

The wall- to bulk-temperature difference divided by the fully developed difference is 
plotted in figure 6 for flow in an insulated channel and parametric values of k. The ratio 
is given by equations (41) and (55) with R set equal to zero therein. The information 
given on this plot, used in conjunction with figure 3, permits evaluation of the wall- to 
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Wall temperature ratio, (t w - t b )/(t 







Dimensionless longitudinal distance, (x/a)/Pe Dimensionless longitudinal distance, (x/a)/Pe 


Figure 6. - Wall temperature ratios for 
internal heat generation in channel 
with insulated walls. Power-law 
velxity exponent, 1/7. 


Figure 7. - Wall temperature ratios for inter- 
na! heat generation in channel with wall 
heat transfer. Power-law velxity exponent, 
1/7. 


bulk-temperature difference at locations along the thermally insulated walls. Thermal 
entrance length is commonly defined as the heated length required for t w - t^ to ap- 
proach to within 5 percent of the fully developed value. A dashed line is drawn in fig- 
ure 6 to facilitate finding the entrance length. For example, for Re = 15 000 and 
Pr = 0.01 (Pe = 150 and k = 1.25 for ^ = 1), the dimensionless entrance length x/a 
is approximately 8, while for Re = 175 000 and the same Prandtl number, Pr = 0. 01 
(Pe = 1750 and k = 2. 50 for ^ = 1), the entrance length x/a is approximately 37. 
Therefore the entrance length increases as the Reynolds number increases. 

The thermal entrance length is also influenced by heat addition at the channel walls 
for various values of k. This is shown in figure 7, which gives the wall temperature 
ratio in the thermal entrance region as a function of parametric values of k and R. 
Then, for example, for Re = 175 000 and Pr = 0. 01, the entrance length, for R = 0. 1, 
is reduced to approximately 32. An approximately 13-percent reduction occurs in the 
entrance length with the given conditions when R varies from 0 to 0. 1. 

In order to assess the accuracy of the present results, it is desirable to compare 
them with existing numerical and experimental data. To my knowledge, no heat-transfer 
measurements for turbulent liquid-metal flows between parallel plates with heat sources 
in the fluid stream are available for comparison. Reference 18 presents an analytical 
study of turbulent heat transfer in fully developed flow between parallel plates with inter- 
nal heat sources. The fluid had a range of Prandtl number from 1 to 100. These results, 
however, are not applicable to liquid metals whose high thermal conductivities give a 
range of Prandtl numbers from approximately 0. 001 to 0. 1. 
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As mentioned previously, it is customary to represent heat-transfer results in terms 
of a Nusselt number. Nusselt numbers in the absence of internal heat generation Nuq 
( given by eq. (44) for k > 1 and by eq. (57) for k = 1) have been evaluated numerically 
for a power-law velocity exponent m of 1/7 and are plotted as solid lines in figure 8 for 
parametric values of k. Also shown in the figure as dashed lines are Nusselt- number 
results for a power- law velocity exponent of zero (slug flow). It is noted that the Nusselt- 
number results for the one-seventh power-law velocity expression are smaller than those 
for the uniform velocity expression. This finding is in accordance with the results for 
parallel plates with uniform wall temperatures and no internal heat sources (ref. 9). 

It is of practical interest to examine the Nusselt numbers in the fully developed re- 
gion, Nuq ( j. These results, given by equation (45) for k > 1 and in equation (58) for 
k = 1 are illustrated in figure 9, where the fully developed Nusselt number is plotted as 
a function of the parameter k. Results for the uniform velocity profile and analytical 
data from references 19 and 20 also appear in figure 9. 

The numerical data of references 19 and 20 are discussed in some detail in refer- 

4 

ence 2. The data in reference 19 is based on experiments in the range 1. 25x10 < Re < 

1. 35x10^, Pr = 0. 0074, while that in reference 20 is given for the range 2000 < Re < 10^, 
0 < Pr < 1.0. The Reynolds and Prandtl numbers to which the Nusselt numbers in these 
investigations correspond were converted to the equivalent values of the parameter k 
(calculated for \p = i// = 1). It is evident that the fully developed Nusselt numbers based 
on the simplified velocity and eddy diffusivity profiles are in good agreement with other 
numerical analyses. 

In view of the foregoing comparison, it may be concluded that the wall-temperature 
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Figure 8. - Nusselt number results for 
uniform wall heat flux without inter- 
nal heat sources. 
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Figure 9 - Comparison of analytical solution 
and numerical data for fully developed 
Nusselt numbers with wall heat transfer 
and no internal heat generation. Average 
eddy diffusivity ratio, 1. 
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results based on these profiles have an accuracy adequate for design purposes and can be 
used to estimate the combined effects of internal heat generation and wall heat transfer. 

It is of interest to consider the effect of internal heat generation on the Nusselt num- 
ber. Consideration will be given here to the fully developed heat- transfer situation, 
thereby eliminating x as an additional variable. If the velocity profile is assumed to be 
uniform, as, for example, at high Reynolds number, then the Nusselt number is inde- 
pendent of the strength of the internal heat source, and the result Nu = Nuq (with Nuq 
given by eqs. (44) and (57)) is valid for all values of the inverse heat-flux ratio S. 

If, on the other hand, the velocity profile is assumed given by a 1/7 power law, then 
the effect on the Nusselt number of internal heat generation is determined from equa- 
tions (42) and (56). Using the numerical data in tables I to IV, the Nusselt numbers 
Nu/Nuq have been plotted in figure 10 as a function of positive values of the inverse heat- 
flux ratio S for parametric values of the diffusivity parameter k. It is evident that the 
effect of internal heat generation is always to decrease the Nusselt number below its 
value in the absence of internal heat generation for all values of k. For given heat- 
source strength, or inverse heat-flux ratio S, the reduction in Nusselt number is 
greatest when transverse heat transfer occurs solely by molecular conduction. The 
presence of turbulent eddying in the liquid metal causes a smaller reduction in Nusselt 
number for a given positive value of S. 



Figure 10. - Variation of fully developed Nusselt number 
with internal heat source strength for 1/7-power-law 
velocity distribution. 
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CONCLUSIONS 


Solutions have been obtained for heat transfer to a liquid metal with an internal heat 
source and flowing turbulently between parallel plates. The power-law velocity model for 
the turbulent velocity distribution was utilized for the analyses; a value of the exponent of 
1/7 was used in the evaluation of the solutions. The effect of transverse eddy diffusivity 
variations has been included through the use of an idealized eddy diffusivity function. The 
wall temperatures and Nusselt numbers in both the entrance and fully developed regions 
can be obtained including the effects of internal heat generation in the fluid and heat trans- 
fer at the channel walls. 

Some of the characteristics of turbulent liquid-metal heat transfer with internal heat 
sources and wall heat transfer and with negligible axial conduction can be summarized as 
follows: 

1. The wall- to bulk-temperature difference is affected by the presence of internal 
heat generation when the velocity profile is approximated by a 1/7-power law. The tem- 
perature difference is unaffected by the presence of internal heat generation, however, 
when the turbulent velocity profile is approximated by a uniform distribution over the 
channel cross-section. 

2. The thermal entrance length for uniform internal heat generation in an insulated 
parallel- plate channel and for a given Prandtl number is increased as the Reynolds num- 
ber is increased. For a given Reynolds number, Prandtl number, and heat-source 
strength, the effect of heat addition at the channel walls is to diminish the entrance length. 

3. The effect of internal heat generation is to decrease the Nusselt number below its 
value in the absence of internal heat generation when the velocity profile is approximated 
by a 1/7-power law. For a given heat source strength and wall heat flux to the fluid, the 
Nusselt number reduction is greatest when transverse heat transfer is entirely by molec- 
ular conduction. 

A comparison of a few of these results with those of others has shown good agree- 
ment and, hence, supports the physical model employed. The results should be adequate 
and useful for preliminary design purposes. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, September 28, 1967, 

129-01-11-05-22. 
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APPENDIX - SYMBOLS 


A hydrodynamic parameter 

A n coefficients in series expansion of temperature for case of internal heat gen- 

eration in insulated channel 

A coefficients defined by product A R (1) 

n n n' ' 

a half-height of channel 

2 2 / 

a Q constant, (k - 1) /4 

B n coefficients in series expansion of temperature for case of wall heat trans- 

fer and no internal heat generation 

B coefficients defined by product B R (1) 

n r n n' ' 

C n coefficients in series expansion of temperature for case of internal heat gen- 

eration in insulated channel 

C coefficients defined by product C Y (1) 

Cp specific heat at constant pressure 

D H hydraulic diameter of channel, 4a 

D n coefficients in series expansion of temperature for case of wall heat trans- 

fer and no internal heat generation 

D coefficients defined by product D Y (1) 

n J r n n' ' 

Ej, E 2 ,Eg arbritrary constants 

F(r) transverse temperature distribution in fully developed region for case of in- 

ternal heat generation in insulated channel 

F(y') transverse temperature distribution in fully developed region for case of in- 

ternal heat generation in insulated channel for limiting value of k = 1 

G(r) transverse temperature distribution in fully developed region for case of 

wall heat transfer and no internal heat generation 

G(y’) transverse temperature distribution in fully developed region for case of 

wall heat transfer and no internal heat generation for limiting value of 
k = 1 


h 

Ij( r ) 


heat transfer coefficient, q /(t - t^) 

integral defined by eq. (27b) 
integral defined by eq. (27c) 
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h<u 

integral defined by eq. 

(27d) 

I 4 (k) 

integral defined by eq. 

(27e) 

I x (t) 

integral defined by eq. 

(50d) 

V 1 ) 

integral defined by eq. 

(50b) 

V 1 ) 

integral defined by eq. 

(50c) 

J l/(m+2) 

Bessel function of the first kind and order l/(m + 2) 

k 

diffusivity parameter, 

Vl + 0. 01 4/ PrRe 0,9 


m exponent in power-velocity expression 

Nu Nusselt number, hDj j//c 

Pe Peclet number, RePr 

Pr Prandtl number, v/a 

Q heat generation rate/volume 

q w prescribed wall heat flux 

R heat-flux ratio, q w /Qa 

Re Reynolds number, 4ua/v 

R(r) function of r 

R n (r) eigenfunction 

r Q radius of pipe 


r^ variable, 1 + (k^ - 1)(1 - y’); 1 < r < k 

S inverse heat-flux ratio, Qa/q w 

o 

s variable, r - 1 

t fluid temperature 

U dimensionless velocity, u/u 

u fluid velocity in x-direction 

u fluid mean velocity 

x longitudinal coordinate measured from channel entrance 

x' dimensionless longitudinal coordinate, 4(x/a)/Pe 

y transverse coordinate measured from channel centerline 



y' dimensionless transverse coordinate, y/a 

z transverse coordinate measured from channel wall, a - y 

a molecular diffusivity for heat, /c/pc p ; dummy variable of integration 

eigenvalue 

/3 dummy variable of integration 

2 

fi eigenvalue 

ejj eddy diffusivity of heat 

eddy diffusivity of momentum 
77 variable, 1 - y’ 

k fluid thermal conductivity 

X 2 n eigenvalue, 4(m + l)a^/(k 2 - l ) m+2 

Arp dimensionless diffusivity, Pr(e^/v) 

v fluid kinematic viscosity 

£ dummy variable of integration 

p fluid density 

r dummy variable of integration 

cp dummy variable of integration 

ratio of eddy diffusivity for heat transfer to that for momentum transfer, / e M 

^ average value of 4/ 

2 9 

eigenvalue, (m + l)/3 n 

Subscripts: 

b bulk condition 

d developed region 

e entrance region 

i entrance value 

Q internal heat generation, insulated wall 
q wall heat transfer, no heat generation 
w wall 

0 either no heat transfer or no heat generation 
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tion because of preliminary data, security classification, or other reasons. 

CONTRACTOR REPORTS: Scientific and technical information generated 

under a NAS A contract or grant and considered an important contribution to 
existing' knowledge. 

TECHNICAL TRANSLATIONS: Information published in a foreign 
language considered to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information derived from or of value to NASA 
activities. Publications include conference proceedings, monographs, data 
compilations, handbooks, sourcebooks, and special bibliographies. 

‘TECHNOLOGY UTILIZATION PUBLICATIONS: Information on tech- 
nology used by NASA that may be of particular interest in commercial and other 
non-aerospace applications. Publications include Tech Briefs, Technology 
Utilization Reports and Notes, and Technology Surveys. 


Details on the availability of these publications may be obtained from: 


SCIENTIFIC AND TECHNICAL INFORMATION DIVISION 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 



